Reduced density-matrix functional theory in quantum Hall systems 
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We apply reduced density-matrix functional theory to the parabolically confined quantum Hall 
droplet in the spin-frozen strong magnetic field regime. One-body reduced density matrix functional 
method performs remarkably well in obtaining ground states, energies, and observables derivable 
from the one-body reduced density matrix for a wide range of system sizes. At the strongly correlated 
regime, the results go well beyond what can be obtained with the density functional theory. However, 
some of the detailed properties of the system, such as the edge Green's function, are not produced 
correctly unless we use the much heavier two-body reduced density matrix method. 



I. INTRODUCTION 

Quantum Hall fluids occur at low temperature in clean 
two-dimensional electron systems exposed to a perpen- 
dicular magnetic fieldi^ Different phases are character- 
ized by the number v that tells the ratio of the num- 
ber of electrons to the number of single-particle states 
in a highly degenerate Landau level (the number of flux 
quanta piercing the sample). Near certain fractional fill- 
ing factors v, the interactions between electrons induce 
an energy gap and lead to a ground state with topologi- 
cal order manifest in exotic quasiparticles and non-Fermi 
liquid edge modest- 
Straight from the outset, numerical simulations have 
been an indispensable guide in development of the 
theory^ While majority of the numerics are exact diag- 
onalization studies only viable with small electron num- 
bers, variational Monte Carlc-^ and density functional 
theor y 13 ' 14 (DFT) have also been applied to larger sys- 
tems. For example, accurate wave functions incorporat- 
ing the complex non-perturbativc effects of electron inter- 
actions have been theoretically devised and later backed 
up by the high overlap with exact numerical results for 
small systems. Due to such developements, reason be- 
hind the energy gap of simplest of the many fractions is 
now well understood within the framework of composite- 
fermion theory that allows for explicit construction of the 
many-body wave function and calculation of topologi- 
cal quantum numbers^ However, in going beoynd the 
compositc-fcrmion theory to more complex phases, the 
Monte Carlo method is crippled since a trustworthy trial 
wave function for the phases we would be interested to 
know more about is not known. In addition, owing to the 
strong correlations, the DFT is inaccurate at, for exam- 
ple, the paradigm fractional quantum Hall state at filling 
fraction v = 1/3 where the vortices supposed to form 
a bound state with the electrons localize at fixed posi- 
tions instead^ Consequently, the exact diagonalization 
is frequently the only viable alternative, leaving the large 
electron numbers beyond the reach of direct calculations, 
though the density matrix renormalization group method 
has brought some progress making a bit larger systems 
computationally feasible ! 16 ' 17 

During the past 10 years, reduced density matrix func- 
tional theory has been revived and applied successfully 



in the chemical physics community* 1 ^ - — The method is 
known to handle e.g. molecular dissociation limits better 
than standard DFT^^ and it has recently been applied 
to homogeneous electron gas. 23 ' 24 In this manuscript, 
we aim to probe the potential of the one-body reduced 
density-matrix functional theory (1-RDMFT) in a frac- 
tional quantum Hall system, specifically a parabolically 
confined quantum dot in the spin-frozen strong magnetic 
field regime. In contrast to many molecular and atomic 
systems where the dominant occupation numbers are typ- 
ically close to one, this is a highly demanding application 
for any numerical method as fractional quantum Hall 
states have long-range quantum entanglement with all 
the occupation numbers small for example near 1/3 in 
the v = 1/3 state. 

The performance of various functionals in predicting 
ground states, energies, and observables attainable by 
the reduced density matrix is compared for small system 
to the exact diagonalization and Hartree-Fock with and 
without the Brillouin-Wigner perturbation theory. For 
larger systems with tens of particles, the comparison is 
done utilizing the accurate Laughlin wave function 4 - for 
filling fraction v = 1/3 state and Monte Carlo methods. 
Energies are produced quite well in all systems at the 
strong-correlation regime v « 1. Even the bulk densi- 
ties appear reasonably good and reproduce the predicted 
edge stripe phase.— However, we are still dealing with an 
approximate method that has its limitations. Detailed 
properties of the edge of the electron droplet, such as the 
edge tunneling exponent ; 26 ' 27 are not produced correctly 
by the present functionals. 

For this reason, we also perform a small comparison 
with the heavier two-body reduced density matrix func- 
tional theory 2 ^ (2-RDMFT) (see the closely related work 
in Ref. HH). Including the exact electron interaction by 
the two-body reduced density-matrix functional appears 
to facilitate a more accurate description of the edge, how- 
ever with a computational cost in large systems beyond 
the reach of present day computers. 

The rest of the manuscript is organized as follows. In 
the next section, we briefly introduce the idea of reduced 
density matrix functional methods. In Section IIII1 we 
describe the model system and derive an exact formula 
for the energy contribution due to one-body operators 
present in our Hamiltonian such that only the interac- 
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tion energy remains to be solved. Section IIV Al and B 
present the details of our 1-RDMFT and 2-RDMFT im- 
plementations, respectively. The 1-RDMFT results are 
divided according to the small or large size of the system 
into Scction lVAland B , followed by the 2-RDMFT calcu- 
lation in Section IV CI Conclusions and future prospects 
arc found in Section [VII 



then reads 



e 

"27 



drdr^-^/t,,,-) 



r - r 



drdr 



,#(r)^(r%(r)&(r' 



(5) 



II. REDUCED DENSITY MATRIX 
FUNCTIONAL THEORIES 

The 1-RDMFT is based on the Gilbert's theorem, 
which guarantees that the ground-state expectation value 
of any observable is a unique functional of the 1-RDM 
7.— It follows that the ground state energy can be writ- 
ten as 



F[ 7 ] = f drdr'<5(r - r')(T(r) + [/(r)) 7 (r, r') + V cc [ 7 ] 

JWl 2d 

(1) 

where T and U are the standard operators for the kinetic 
energy and an external potential while the functional for 
the interaction energy 14 [7] is unknown. It is simple to 
show that this functional yields the exact ground state 
energy for the exact 1-RDM 



7 (r,r') =7V J **(r,r 2! r 3 ,...,rjv) 



(2) 



x *(r', r 2 , r 3 , . . . , r/v)dr 2 dr 3 . . . dr 



JV 



The form of / could vary in different type of systems, 
whereas those used in this study are enlisted in Table U in 
Sec. IIV Al In analogy with the density functional theory, 
the first term is referred to as the Hartree term while the 
latter is the exchange-correlation term. However, com- 
pared to the DFT, 1-RDMFT has a couple of advantages. 
Firstly, the method obtains not only the density but the 
whole 1-RDM so that the ground-state expectation value 
of any one-body observable can be readily computed. 
Thus, for example kinetic and interaction energies can 
be readily separated and Green's functions calculated. 
Secondly, although both methods are in principle exact 
for an exact functional, due to the variable 7 (vs. p), 
it is easier to develop a good 1-RDM functional than a 
good density functional. This is the reason why the DFT 
does not work in the strongly correlated regime where the 
proper treatment of electronic correlations is important. 
However, there are fresh ideas of how to treat strongly 
correlated electrons with DFT i 32 ' 33 



III. QUANTUM DOT MODEL 



if V ee is replaced by half the Coulomb energy of the exact 
pair-density 



E ee 4f drdr'^4 
2e . B 2d r - r 



p 2 (r,v')=N(N-l) J **(r,r',r 3 ,...,r w ) 
x *(r, r', r 3 , . . . , r N )dr 3 . . . dr N 



(3) 



The 2-RDMFT minimizes the resulting exact functional 
F subject to a subset of complete W-representability con- 
ditions, known as the 2-representability conditions, and 
other possible constraints due to additional symmetries 
(see Sec. IIVB[) . On the other hand, the crux of the 1- 
RDMFT is to approximate the pair-density by a func- 
tional of the l-RDMi£L The customary way to do the 
approximation, which we will also employ in this paper, 
is to replace the pair-density above by 

p(r)p(r')-J2f(ni,nM(^A^)^)Mr') (4) 

where p(r) = 7(1", r) is the density at r, <fii are the natural 
orbitals (eigenvectors of 7), and / is a function solely of 
the natural occupation numbers n, € [0, 1] (eigenvalues 
of 7). The functional for the interaction energy in Eq. (UJ) 



The quantum Hall droplet is modeled by the two- 
dimensional effective-mass Hamiltonian 



H = J2 



2m* 



E- 



(6) 



where N is the number of electrons, A is the planar 
vector potential of the homogeneous magnetic field B 
perpendicular to the sample plane, and the energy scale 
of the external confinement potential hioo is typically a 
few me Vs. The material parameters are the effective 
mass of the electrons m* = 0.067m e and the dielectric 
constant of the GaAs semiconductor medium e = 12.7. 
Coulomb interactions tend to spontaneously polarize the 
electron spins in an effect known as quantum Hall ferro- 
magnetism. For this reason, the relatively weak Zeeman 
term has been omitted and the electrons are assumed 
spin-polarized. From here on, we use oscillator units so 
that the energies arc expressed in units of hw and lengths 
in units of \Jh/m*ui where lo 2 — loq + (u c /2) 2 with cy- 
clotron frequency w c = eB/m*c. 

For N = 1, the energy states are written as 



7r(n + to) 



■z m m{zz)e 



-zz/2 



m > —n, n > 



(7) 
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where z = x + iy and are the generalized Lagucrre 
polynomials. The corresponding eigenvalues are given by 



E': 



2n+l 



in 



(8) 



We take on interest in developing a computational 
method for the fractional quantum Hall states at the 
strong magnetic field regime, so we may assume that 
ujq <C oj c . Then Sw c /2 is close to unity, such that the 
term in parentheses in Eq. ((SJ) becomes small and values 
of quantum number n other than become irrelevant to 
the low energy physics. This is the Landau level projec- 
tion to the band with n = 0. It should not be difficult 
to generalize the 1-RDMFT method to include spin and 
several Landau levels and study the region v > 1 as well, 
however, for simplicity we stick to the spin-polarized low- 
est Landau level v ^ 1 in this study. 

Note that in our two dimensional model, m in Eq. ([7]) 
is the one and only angular momentum quantum number. 
We can write the contribution to the total energy due to 
terms other than the Coulomb interaction for a system 
of N electrons in the lowest Landau level n = in terms 
of the total angular momentum (quantum number) M 
exactly as 



T + U 



N r 

i=l L 



1 



m, = N+ 1 M 



~2 



(9) 

As the total angular momentum operator M commutes 
with the total Coulomb interaction energy operator V ee , 
to solve the Landau level projected spectrum, the remain- 
ing task is to find the common eigenstates of the Landau 
level projected V cc and M . 

For TV and M small enough, these are solved exactly 
with the configuration interaction method (exact diag- 
onalization) since the number of possible single-particle 
states is finite. We may go to a bit larger N and M 
by taking interest in only the ground state and finding 
it by the Lanczos algorithm (cf. Refs. |36||37| ). However, 
the exponential growth of the many-body basis limits the 
particle number to around 10 depending on M, and some 
kind of an approximate calculation method becomes nec- 
essary. The Kohn-Sham DFT is still a good method for 
the weakly correlated regime^ but for larger M the nat- 
ural occupations tend far from 1 and the DFT no longer 
gives us good results. Absence of a reliable trial wave 
function for generic M makes the implementation of vari- 
ational Monte Carlo rather uncertain ! 12 ! 35 In this paper 
we are going to see, if the reduced density matrix func- 
tional theory can make itself useful. The expectation is 
that it will work better than DFT at least when the eigen- 
values of the density matrix are small meaning i/ « 1. 
With several Landau levels, this would correspond to the 
situation where the highest occupied Landau level has 
low filling fraction. 



IV. COMPUTATIONAL METHODS 
A. 1-RDMFT 

Due to the exact formula for the one-body operators' 
energy contribution (Eq. ©), the energy functional of 
Eq. (Q} simplifies to 



F[y] = N+[1 



M + V ee [-y] , 



(10) 



with the constraints that the particle number and the to- 
tal angular momentum are N and M , respectively. More- 
over, the natural orbitals for energy state |^) in the low- 
est Landau level are directly the single-particle energy 
states since the angular momentum conservation yields a 
diagonal density matrix in this basis 



In 



(ii) 



where a m and its adjoint annihilate and create a particle 
with quantum numbers n — and m (see Eq. IJT])). As 
N — 1 particles have at least a total angular momentum 
(N — 1)(N — 2)/2, the maximum single-particle angular 
momentum becomes k = M — (N — 1)(N - 2)/2. There- 
fore, the minimization of the interaction energy reduces 
to the constrained minimization of a function whose vari- 
ables arc k + 1 occupation numbers 



KetMLo) 



2 ^ 



(runjVijij - f(ni,nj)Vijji) (12) 



where V^-jy are the interaction matrix elements of the 
lowest Landau level orbitals computed in Ref. [HI, and 
the constraints are explicitly written as 



k k 

n-m = N and 

m— m— 



M 



(13) 



To optimize the occupation numbers, we first express 
them in terms of variables 9i such that rij = sin 2 6i . 
The minimization of V cc with the above two remaining 
constraints is then performed with the interior point or 
Neldcr-Mead algorithm of Mathematical 

A vast number of functions /(n^n^) have been pro- 
posed in the literature in treatment of simple atoms and 
molecules, and it is not immediately clear, which of them 
would work well in the current system. Table U summa- 
rizes those used in this work for small systems to find 
out the optimal ones. For large systems, we only use 
the density- matrix power functional /(n;, rij) = (rii, rij) a 
(P-a)i^i Each form of the off-diagonal f(rii,njjti) corre- 
sponds to two different functions, first of which has diag- 
onal f{fii,riij = rii while the second has f(ni,rii) = ?i 2 . 
Since the natural orbitals are orthonormal, integration 
of the approximate pair-density in Eq. (£[]) over the co- 
ordinates yields the proper normalization N(N — 1) for 
f (rii, rii) — rii. The form n 2 , on the other hand, is jus- 
tified as it cancels the self-interaction terms Van arising 
from the Hartrec term. 
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A few words about different forms of the off-diagonal 
terms are in order (see last column in Table H]) . The form 
nfn^~ a for / was first introduced by Miiller, who found 
a = 1/2 (MU) to be the optimal value.— Goedecker and 
Umrigar (GU) used the modification to the diagonal that 
removes the self-interaction terras^ Much similar to the 
density-matrix power functional, we generalize these to 
MU-a and GU-a where the square root is replaced by an 
arbitrary power, presumably lying between half and one. 
Alternatively, it is often physically motivated to reduce 
the overcorrelation of MU by switching the sign of the off- 
diagonal terms between weakly occupied natural orbitals 
W (BBCl) or additionally banishing the square root for 
pairs of strongly occupied orbitals S (BBC2)i 22 i 23 For 
simplicity, we define the strongly and weakly occupied 
orbitals directly from the Hartree-Fock solution where 
the occupations are or 1. 



TABLE I: Functions f(rn,nj). S and W refer to the strongly 
and weakly occupied natural orbitals, respectively. 





f(rn,Tu) 




MU 
GU 


rii 
n? 




MU-a 


m 






GU-a 


ni 




(mnj) a 


P-a 


nt a 




BBCl 
BBCIS 


m 
n\ 




I -^/niUj i,j G W 
1 l/rnnj i,j iW 


BBC2 
BBC2S 


rii 
nf 




— y/rnrij i, j £ W 
riiUj i,j e S 
yfrnrij else 



Finally, we would like to point out an issue with the 
physicality of the obtained solution, which is generally 
more of a problem in the higher order RDM methods. 
Specifically, Coleman has shown that necessary and suf- 
ficient condition for the 1-RDM to be 7V-representablc, 
meaning that there exists 1^) in the Hilbert space of the 
system such that 7^ = {^\a\aj\^), is that its eigenval- 
ues are between and 1 and their sum is A r 4i However, 
if we additionally demand a symmetry, it may be that 
the obtained solution is not rcpresentable in the symme- 
try restricted part of the Hilbert space. To be explicit, 
if in our model we demand total angular momentum of 
a two electron system to be 2, the symmetry restricted 
Hilbert space of states with m = 0, 1, and 2 consists 
only of one state with occupations (1,0,1) while 1-RDM 
method could give us unphysical occupations (0.5, 1, 0.5) 
both having the same particle number and angular mo- 
mentum (0.5 x0 + lxl + 0.5x2 = 2). While this could 
be avoided by imposing additional constraints, we refrain 
from doing it as that would be exponentially unfeasible 
for larger systems. Additionally, it is plausible that the 
unphysical solutions have less weight when the size of the 
physical Hilbert space increases. 



B. 2-RDMFT 

In analogy with the 1-RDMFT, we now minimize a 
functional of the 2-RDM 

T% = (tt|otota»o fc |*> . (14) 

The marked difference is that we now have an exact func- 
tional for the interaction energy 

^ee(r) = \ ^Jkl , (15) 

i,j,k,l 

however, with the cost of large number of additional pa- 
rameters to optimize with similarly large number of ad- 
ditional constraints. Furthermore, the set of constraints 
to be listed below form a relatively stringent set of nec- 
essary conditions that only in special cases is sufficient 
for the obtained solution to be exactly ./V-representable, 
meaning a physical l^) to exist such that Eq. (fl"4|) holds. 

The minimization of the functional (|15[) is performed 
by forming an augmented Lagrangian function and min- 
imizing it following the algorithm in Ref. |42| . The 
minimization is performed with limited memory quasi- 
Newton algorithm of Mathematica, and the form of the 
augmented Lagrangian is 

L = F[ 1 ]-J2^c t + J2cl/^ (16) 

i i 

where are the Lagrange multipliers and \i > is the 
augmentation parameter used to enforce the convergence 
of the constraints c,; = 0. In the following, we first intro- 
duce the subset of applied TV-representability conditions, 
and after that, impose the further constraints due to the 
fixed total angular momentum, M-representability. 

1. N -representability 
The trace condition 

2>:i = ffi_^ (17, 

i<j 

is used to fix the particle number to N. 

Positivity conditions form the major part of the N- 
representability constraints. Consider an operator of 
the form A = V, (j , > <•■ • • ■ Since 
> 0, it follows that 

For k = 2 we obtain the 2-positivity condition for the 
2-RDM 

ij 
hi 
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By a different choice of A, positivity conditions of the 
exact same form can be derived for the two other repre- 
sentations of the 2-RDM 



Qki =<*l«i<»J*i4l*> and 
G% =<tt|oto,-aJo*|tt) . 



(20) 



While the representations T, Q, and G are all equiva- 
lent as they are related by the fermionic anticommutation 
rules, the positivity conditions are incquivalent and must 
be taken into account simultaneously. We use the anti- 
symmetric basis = — \ ji))/2 for T and Q ma- 
trices since T% = -T J k \ = -rS and Q« = -QJJ = -QH 
Furthermore following Ref. |42j, since all the three matri- 
ces are real and symmetric under ij o kl, the positive- 
definite condition can be accounted for simply by writing 
the matrices as square of symmetric matrices T = R 2 , 
Q = S 2 , and G = T 2 (meaning rgj = E( M ) 
etc.) and optimizing the upper diagonal of matrices R, 
S, and T. The linear relations linking T and Q and T 
and G become the relevant constraint equations 



Q 13 



Sjkjii + SikSji - SuSjk , 



(21) 



where the 1-RDM 7^ may be obtained through 

1] Gwb = Nlij or £ ^ = (N - lfry . (22) 

2. M -representability 

Since m i a l ct il 1 ^ r ) = -^l 1 ^)) we can form one inde- 
pendent nontrivial equation involving 2-RDM (equivalent 
to the contracted Schrodinger equation in Ref. 1291 ) 



^m k G 



(23) 



The trace of this is already fixed by Eqs. (|17l) and (|2"21 



(24) 



Since in our system the angular momentum quantum 



numbers 



= j are always non-negative, the maxi- 



mum angular momentum for a pair of electrons to have is 
k 2 = M - (N - 2)(N - 3)/2 where the subtracted term is 
the minimum angular momentum of TV — 2 electrons. Ad- 
ditionally, some of the matrix elements of the 2-RDM are 
zero because states with different total angular momen- 
tum are orthogonal. The independent constraints due to 
these considerations read 



1 kl 







■ j ^ k + 1 or i + j > k% 



(25) 



while Eq. (|2T[) communicates them to Q and G. Though 
we can just drop the corresponding terms from our equa- 
tions, we still need to take into account the ensuing con- 
straints on the actual variables R.S. and T. 



C. Monte Carlo 

The calculation of the natural orbital occupations is 
relatively simple using the Monte Carlo technique. Un- 
like in the previous Monte Carlo study in Ref. HH, we 
know from the start the natural orbitals that diagonalize 
the density matrix, and thus we only need to calculate 
the occupations. 

Starting from the definition ©, the orthogonality of 
the natural orbitals, and the expansion of 1-RDM using 
the natural orbitals 



7 (r,r') = Y^n m (f>* m {r)(j) m {r') , 



the occupations are integrated as 

n m =N J **(n,r 2) ...)*(r',r2 
x (^ m (ri)</4(r')dr'tfri 



■ •) 



dr 



N 



This is further reformulated as 

n m =N 



J l*(ri,r 2 , 



*(r',r 2) ...) 



*(ri,r 2 ,.. 
x m (ri)#£,(r')dr'dri . ..dr N 



(26) 



(27) 



(28) 



which can be symmetrized and rewritten in Monte Carlo 
expectation value as 




*(r',r 2 ,.. 



-<M ri )C(r')dr' 



*(ri,r a ,...) 

' {ri}e|*| 2 

where the summation is over the coordinates {r^}^. 
The first strategy for Monte Carlo evaluation of the oc- 
cupation is to sample these coordinates from |\&| 2 and to 
integrate over r' on a grid 42 

For a better option in our case, we first rewrite 

m (r*)0m(r') 

Carlo integration also over r' as 



I0m(r')| 2 0™(ri) and then do a Monte 



E 



*(r',r 2 , 



0m (r,) 



*(n,r 2 ,...) m (r') 



(29) 

where {rjj-^Lj is again sampled from |^| 2 and r' from 
\4> m \ 2 - This option can be made more stable by noting 
that the natural orbitals have rotation symmetry and 
1 0m (r') 1 2 depends only on the radial coordinate r' and 
not on the angle 9' . Now, the radial integral over r' can 
be done using Monte Carlo integration and the angular 



integral by averaging over a uniform grid 
*(r',r 2 ,...) & n {n,6 t )\ 



as 



1 X - * 

hi 



(r l!r2! ...)0m(r',^) 



{r;}e|*| 2 , r'e|0„ 
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Notice that r' is generated separately for each <p. 



V. RESULTS 

In the following, we will first compare the performance 
of various 1-RDM functionals in calculating the interac- 
tion energies in different (N, M) sectors in exactly solv- 
able small quantum Hall droplets. This analysis is deep- 
ened by analyzing the occupation numbers obtained for 
the ground states. Once we have established that we have 
a decent functional, we proceed to test its performance 
in larger systems. Since no exact results are available for 
larger systems at the strongly correlated regime i/« 1, 
we employ the next best gripping handle, which is the 
Laughlin's variational wave function for filling fraction 
v = 1/3. Moreover, we compare the energies and occu- 
pations obtained with density-matrix power functionals 
(P-a) to the results extracted from the Laughlin's wave 
function with Monte Carlo techniques. From the occupa- 
tion numbers we calculate the edge Green's function Ge, 
which gives information about the topological order of 
different quantum Hall phases. It is interesting to see if 
the newly applied methods (1-RDMFT and 2-RDMFT) 
reproduce the correct decay properties of Ge- 



A. Towards a good 1-RDM functional in small 
quantum Hall droplets 

The energy states of the quantum dot may be written 

as 

E n MN = Nhw + (Hu- hw c /2)M + V c n c (M, N) , (30) 

where y™ (M, N) is the nth eigenvalue of the total inter- 
action energy in sector (M, N), whose dependence on the 
physical parameters is scaling by a factor e 2 /el. Because 
of the second term, it follows that all possible ground 
states for N particles are found at the intersections of 
the M-V ee -c\xrve and its convex lower envelope. 

Figure QJa) shows the M-V^-curves for 4 elec- 
trons computed with the configuration interaction (CI), 
Hartree-Fock method (HF), and the Brillouin-Wigner 
2nd order perturbation theory (BW) to the HF state. 
The exact diagonalization (CI) ground states, detected 
by the convex envelope, are marked by circles. While 
both HF and BW predict the correct ground states, the 
perturbation theory leads to a significant improvement 
to the HF energy. The 2nd order perturbation theory 
is very accurate and close to the CI result for M < 14, 
and it is roughly half-way between HF and CI energy for 
M > 14. 

Physically the cusp structure seen in the results fol- 
lows from the energetic advantage of a configuration, 
in which the 4 electrons are located at vertices of a 
square. This configuration has non-zero weight only if 
the angular momentum attains a special value such that 



N(N - l)/2 = mod(M,N). The difference between sub- 
sequent magic angular momentum states is the number 
of vortices found at the center of the system. As mag- 
netic field is increased, vortices that carry quantized an- 
gular momentum, and in a sense quantized magnetic flux, 
emerge at the center. As the particle number is taken 
very large, the number of cusps increases while all cusps 
no longer correspond to a ground-state at certain sys- 
tem parameters. Instead, a few of them are more special 
than others forming the incompressible vacuum state of 
certain fractional quantum Hall phase as the remaining 
cusps are related to the quasiparticle/vortex excitations 
that are eventually responsible for the finite extension of 
the quantum Hall plateaus. 

The corresponding energies obtained with 1-RDMFT 
are shown in Figs.QJb) and (c). At this point the choice 
of parameters a = 0.75, a = 0.7, and a = 0.65 for MU- 
CH, P-a, and GU-a is an educated guess, whereas the 
effect of the parameter a becomes apparent in the next 
section (basically, it tunes the strength of electron cor- 
relations). Typical to 1-RDMFT calculations in general, 
the energies are below the CI result. The dashed lines 
systematically lie below the solid lines of the same tone, 
due to the self-energy cancellation present in the latter. 
Basic MU and GU functionals clearly overestimate the 
correlation energy and behave even qualitatively wrong 
as they fall too fast with increasing M. For the rest of 
the functionals, V cc seems to decline at about the correct 
rate as a function of M. However, a nice cusp structure 
is only seen with the BBC1S and BBC2S functionals, 
which inherit the cusps from the HF state used in the 
selection of the strongly and weakly occupied orbitals. 
Overall the energetically best of these 1-RDM function- 
als (P-0.7, GU-0.65, BBC1S, and BBC2S) perform better 
than the 2nd order perturbation theory when M > 14, 
although only BBC2S produces the correct ground state 
structure. 

The equivalent curves for 6 electrons are shown in 
Figs. [T] (d)-(f). As the 6-electron results are quantita- 
tively like the 4-electron results, the performance of dif- 
ferent functionals seems to be rather insensitive to the 
particle number. For six electrons, the cusps should oc- 
cur at N(N - l)/2 = mod(M,JV) or N(N - l)/2 = 
mod (M, AT — 1) corresponding to a hexagonal configura- 
tion or a pentagonal configuration with one electron at 
the center. All the cusps are not correctly reproduced 
with any 1-RDM functional, though BBC2S result fol- 
lows the cusp structure quite well at M > 27. 

The quantum Hall droplet model has the property that 
an increase in M increases the area of the droplet and also 
the electronic correlations quantified in reduction of the 
occupation numbers of the natural Landau level orbitals. 
The average occupation number of the relevant orbitals 
is close to the corresponding macroscopic quantum Hall 
filling fraction v and becomes exact as N is taken to 
infinity. For example, the v — 1 state occurs at the 
minimum angular momentum M — N(N — l)/2 where 
the first N orbitals have occupation 1. Second example is 
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FIG. 1: (Color online) The minimum interaction energy V ee at each angular momentum M for 4 (a-c) and 6 electrons (d-f). 
The exact diagonalization ground states detected by the convex envelope are marked with circles in (a) and (d). 
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FIG. 2: (Color online) The occupation numbers for selection of the 4-electron (upper panel) and 6-electron (lower panel) ground 
states indicated in Fig. Q] (a) and (d), respectively. A bar that fills the assigned height corresponds to occupation number 1. 
Occupations are ordered according to the increasing single-particle angular momentum m of the lowest Landau level orbitals 
starting with at the left. Because most likely location of electron at orbital m is at distance r = \pm from the center, one 
can think the set of orbitals as a radially discretized disk. 



8 



the Laughlin's wave function^ for filling fraction v = 1/3 
*L /3 (W) = Il^-^) 3 e-' E ^^ • (31) 

It has angular momentum M = 3N(N-l)/2 and 3(N - 
1) + 1 nonempty orbitals such that on average the fraction 
v = N/(3(N - 1) + 1) ps 1/3 is filled. The Laughlin state 
has about 0.98 overlap with the exact ground state for 4 
and 6 electrons, and it is the lowest angular momentum 
zero-energy state of the short-range model interactional 



(32) 



Note, however that while for four electrons the highly cor- 
related Laughlin state occurs at M = 18 near the center 
of the M-window in the M-V ee - curves, the corresponding 
angular momentum for six electrons is M = 45, and it 
is the highest M included in the corresponding figures. 
Overall, it seems that the perturbation theory works en- 
ergetically well near v = 1 where the correlations are 
weak while the 1-RDM functionals perform better at the 
strongly correlated regime v <C 1, which raises some hope 
for the 1-RDMFT to prove valuable in quantum Hall sys- 
tems. 

But how close are the obtained minimizing 1-RDMs 
actually to the exact results? Recall that the natural 
orbitals in the lowest Landau level are fixed and their oc- 
cupations completely specify the 1-RDM. Figure [2] shows 
the occupation numbers corresponding to the ground 
states indicated in Fig. [1] (a) and (d) . 

Looking at the first row of occupation numbers for 4 
electrons, the next HF ground state is obtained from the 
previous by adding a hole to the center leading to angu- 
lar momentum increase TV. The exact CI result below 
is similar but, in addition, the correlations spread the 
occupations at each step. On the third row, the sec- 
ond order perturbation theory BW has a small spread 
of occupations in accordance with the energy curves in 
Fig. HJa). The 1-RDM functional results on the follow- 
ing 9 rows are varied in nature. In accord with the poor 
energies, MU functional leads to a way too large spread 
of occupation and so does the GU, although the latter 
also pins some occupations to one. The inclination to- 
wards pinning is due to the self-energy cancellation, since 
without the cancellation non-pinned occupations lead to 
negative self-energy contribution lowering the total en- 
ergy (n| — Hi < for < rii < 1). This is the reason why 
MU-0.75 is more spread out than GU-0.65, and BBC oc- 
cupations are a bit less pinned than BBCS occupations. 
Despite the better energetics, the GU-0.65, BBC1S, and 
BBC2s occupations seem not to be much better than the 
second order perturbation theory. On the contrary, P- 
0.7, on the other hand, has both quite good energy and 
occupations numbers only slightly less spread than the 
exact result. Note however, the nonzero first occupation 
at M = 26 and M = 30 in contrast to the exact result. 

For six electrons (Fig. [Tj (lower)), the occupations are 
similar. Note the high probability for one electron to be 



at the center in some of the HF and CI ground states, 
correctly reproduced by many of the functionals. The 
six-electron occupation numbers at M = 30 and M = 35 
can be directly compared to the occupations obtained 
with DFT in Ref. [l5|, and they are found to be a bit 
similar to our BBC1 or BBC2 results. 

On the whole, the P-0.7 power functional seems like a 
good candidate functional for systems with large number 
of electrons. GU-a with a < 0.65 and MU-q with a > 
0.75 could also work well, however, since the diagonal 
part of the P-a functional is somewhat a compromise 
between these two, we employ the P-a functional in the 
remainder of the manuscript. 



B. 1-RDM at large N 

The results of the previous section suggest that the 
density matrix power functional (P-a) could be a good 
functional in quantum Hall systems, and thus we ap- 
ply it to large systems for a few parameters a. More- 
over, we concentrate on the v = 1/3 state, whereby close 
to exact nonperturbative results can be computed with 
the Laughlin's variational wave function (Eq. (|31[1 ) using 
Monte Carlo. It is natural to limit the number of natural 
orbitals to that of the Laughlin wave function 3(N— 1) + 1 
although the realistic Coulomb ground state would actu- 
ally extend, weakly though, to a few more orbitals. 

The occupation numbers obtained in such way for 
-/V = 10, N = 20, and N = 30 are presented in Fig. G2 
As seen in the exact result (first row in Fig. G3a)) the 
long-range Coulomb correlations cause oscillations in the 
occupations around 1/3 (CI), apart from the edge den- 
sity modulation not present in the Laughlin's occupa- 
tion numbers (MC). Sliding a from 0.6 to 0.75 gradu- 
ally strengthens these oscillations. P-0.65 is close to the 
Laughlin's occupations while P-0.675 is close to the exact 
occupations (CI). Similar behavior is seen at larger par- 
ticle numbers in Figs. [3](b) and (c), where P-0.675 yields 
again occupations plausibly closest to the exact unknown 
result. 

The oscillations in the occupation numbers reflect the 
formation of an edge striped phase. An extrapolated 
phenomcnological formula for the slow-decaying charge 

v = 1/2 
lei. M 



density oscillations at the v = 1/3 edge at the thermo 
dynamic limit is given in Ref. 



p( S ) = i(Erf( S ) + l) 



1)) 



(33) 



where s/y/2 is the distance from the edge located at 
y/3(N — 1), Erf is the Gauss error function, and Jo is the 
Bessel function of the first kind. Fig. 0Ja) shows the ra- 
dial charge densities (r)i/j(r)\^) calculated from the 
30 electron occupation numbers compared to the extrap- 
olated formula (red dashed line). The latter has slightly 
longer oscillation wavelength compared to the P-a results 
while the amplitude of oscillations suggests that the opti- 
mal value of a is somewhere between 0.675 and 0.7. Thus, 
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FIG. 3: (Color online) The occupation numbers at the angular momentum of the 1/3 Laughlin state for (a) TV = 10, (b) 
TV = 20, and (c) TV = 30 electrons with 3(TV — 1) + 1 natural orbitals. MC is the exact Laughlin wave function result extracted 
with Monte Carlo. The dashed line marks 1/3 occupation. 
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FIG. 4: (Color online) (a) 30 electron density profiles corresponding to Fig. EJc) shifted vertically by 0.5 unit. Red dashed line 
is the phenomenological estimate for the density oscillations at the thermodynamic limit. (b),(c) The decay of the edge Green's 
function calculated at r — \/3(N — 1) + 1 from the occupations in Fig. [3ja) and (b), respectively. The short red lines illustrate 
slopes —3 (theoretical prediction for v — 1/3) and —1 (Fermi liquid). 



although Eq. ([55)1 has zero free parameters, it matches 
the 1-RDM results reasonably well. 



Edge Green's function Ge is the amplitude for electron 
to propagate a distance along the edge. In the quantum 
Hall droplet, the distance is related to the angle 9 be- 
tween the two points, and Ge = (^Itp^ (zoe l6 )tp(zo)\ , $) 7 
where z is a point of the edge. Chiral Luttinger liquid 
theory of the fractional quantum Hall edge predicts the 



universal asymptotic behavior— ^ 

|G E | cx \z Q e i9 - ZqI-o cx | sin(0/2)r 9 , (34) 

where g = 3for^ = l/3. Values of g 7^ 1 lead to non- 
Ohmic current- voltage dependence / cx V 9 in the tunnel- 
ing experiments. However, thus observed experimental 
value, g as 2.2 — 2.8, is contrary to the theory possibly 
sample dependent ! 26 ' 27 

The decay of |Ge| calculated with the density- matrix 
power functionals is compared to the exact (CI) and the 
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Laughlin wave function's result (MC) in Figs. @|b) and 
(c) for TV = 10 and TV = 20, respectively (not shown case 
TV = 30 looks similar). The short lines signify power-law 
exponents g = 1 and g = 3. Except for the curve cor- 
responding to P-0.75, which oscillates heavily, the P-a 
curves follow closely the theoretical dashed black line un- 
til sin (0/2) « 0.3, after which they sheer off the course to 
yield an exponent 1. This is due to the incorrect weights 
of the occupation numbers near the edge of the system 
and investigated further in the next subsection where we 
apply 2-RDMFT to a smaller system. 

The interaction energies are shown in Table [TTJ The 
exactness of the Laughlin trial wave function's energy 
(MC) up to 0.1% for TV = 10 is expected to carry on to 
the larger electron numbers. The interaction energy is 
seen to increase as a function of a in P-a and is optimal 
with the functional P-0.7, which attains 99.9% of the 
interaction energy for TV = 20 and 30. Typical to 1- 
RDMFT calculations in general, the energies are mostly 
below the assumed nearly exact Monte-Carlo energy. 

While knowledge of the ground state energy may be 
useful when comparing different methods, only energy 
differences are physically meaningful. In the RDM meth- 
ods, we can calculate the energy differences between 
lowest energy states of different (M, TV) sectors such as 
addition energy (change in TV) as well as quasiparticle 
and some edge excitations (M changes). If the excited 
state becomes the ground state for some parameters, the 
Gilbert's theorem guarantees the existence of a 1-RDM 
functional minimized by the exact 1-RDM but even if 
this is not the case, a good functional might still exist. 

As mentioned previously, some of the cusps in M- In- 
curves correspond to the quasiparticle excitations of sta- 



ble quantum Hall phases. Next, we consider such a quasi- 
hole excitation above the v = 1/3 state. The quasiholc is 
a charged vortex carrying fractional charge q = e/3 and 
obeying anyonic statistics To obtain interaction part 
of the quasihole excitation energy at v = 1/3, we need to 
calculate the difference V ee (M 1/3 +N, TV) - V cc (M 1/3 , TV). 
The angular momentum M\ + TV follows from the 
Laughlin's quasihole wave function, which is also used to 
compute an estimate for V ee (M 1 / 3 + TV, TV) with Monte 
Carlo. Due to the fact that Laughlin's wave function is 
more accurate than the quasihole wave function, varia- 
tional principle implies that the Monte Carlo estimate 
to the (negative) contribution to the excitation energy 
is likely an upper bound to the exact result. Neverthe- 
less, for 8 particles the difference to exact CI result is 
less than 0.2% so we expect the estimates to be quite 
accurate. The quasihole wave function reads 

*l /3 ({^» = n^-^n^-^) 3 ^ 1 ^ . (35) 

i i<j 

where zqm = (V-^0 Si z «- The interaction contribution 
to the excitation energy is shown in Table ILTI1 The power 
functionals appear to overestimate the energy gap by 10 
to 20 percent compared to the trial wave function though 
the results seem to get more accurate with increasing TV. 
This preliminary result indicates that the method could 
prove useful in assessing the stability of different models 
for quantum Hall phases characterized by certain angu- 
lar momentum and spin. Additionally, 1-RDM method 
offers a simple framework to include the higher Landau 
levels, however, instead of the bare eigenvalues of the re- 
duced density matrix, one must then optimize the eigen- 
vectors also. 



TABLE II: Interaction energy in units of e 2 /el at angular 
momentum M = 3TV(TV — l)/2 and the percentage captured 
of the energy of the Laughlin wave function for density- matrix 
power functionals P-a. 





TV = 10 


TV = 20 


N = 30 




N = 10 


TV = 20 


TV = 30 


CI 


10.14 








99.9% 






MC 


10.15 


32.92 


64.00 




100.0% 


100.0% 


100.0% 


P-0.6 


8.93 


30.36 


60.08 




88.0% 


92.2% 


93.9% 


P-0.65 


9.58 


31.72 


62.15 




94.4% 


96.4% 


97.1% 


P-0.675 


9.87 


32.34 


63.10 




97.2% 


98.2% 


98.6% 


P-0.7 


10.12 


32.88 


63.93 




99.7% 


99.9% 


99.9% 


P-0.75 


10.44 


33.44 


64.80 




102.9% 


101.6% 


101.2% 



TABLE III: Interaction energy of one elementary quasihole 
excitation in units of e 2 /el and the percentage captured of 
the interaction energy of the wave functional quasihole model 
(MC). 





TV = 10 


TV = 20 


TV = 30 




TV = 10 


TV = 20 


TV = 30 


MC 


-0.314 


-0.497 


-0.635 




100% 


100% 


100% 


P-0.6 


-0.396 


-0.601 


-0.758 




126% 


121% 


119% 


P-0.65 


-0.381 


-0.588 


-0.736 




121% 


118% 


116% 


P-0.675 


-0.376 


-0.583 


-0.733 




120% 


117% 


115% 


P-0.7 


-0.376 


-0.584 


-0.689 




120% 


118% 


109% 


P-0.75 


-0.378 


-0.585 


-0.326 




120% 


118% 


51% 



C. 2-RDMFT results for three electrons 

In this final part, we apply the exact 2-RDM functional 
(Eq. ([3])) to a three electron droplet in the 1/3 state again 
with the maximum single-particle angular momentum set 
to 3(TV - 1). We will see that the resulting 2-RDM, 
though not strictly physical, is close to the exact solution 
and yields better results than our 1-RDM functionals. 

Compared to the 1-RDM calculations seen above, the 
computational cost of the problem in the 2-RDM opti- 
mization is considerably larger and scales at higher order 
p e (versus p 4 ) with the number of single-particle orbitals 
p. The 1/3 Laughlin state for 3 electrons has angular 
momentum M = 9 and requires only 7 single-particle 
states. However, the number of optimization variables in 
L, Q, and G are 276, 378, and 1225, respectively, and the 
constraint equations (|17l21l23l25p together lead to 2798 
constraints each facilitated by a Lagrange multiplier. In 
practice, this means that this method takes more time 
than exact diagonalization in any system that could be 
solved in a reasonable time. However, owing to the ex- 
ponential scaling of the exact diagonalization problem, 
the situation could change in future with development 



11 



of faster computers and more efficient semi-definite pro- 
gramming and optimization algorithms. 

Figure [5] shows the occupation numbers and the decay 
of the edge Green's function for three electrons at v = 1/3 
(AI = 3N(N - l)/2 = 9). Due to a finite size effect, the 
exact diagonalization Green's function has a downward 
cusp at sin(0/2) w 0.9. The 2-RDMFT result has a sim- 
ilar cusp while the P-0.675 result, which turns smoothly 
to a lower slope decay, does not have one. We verified 
that this is due to the difference in the weights of the last 
three occupation numbers corresponding to the edge of 
the system. Since the HF solution would have exponent 
g = 1, the correct decay property of the Green's function 
follows from the off-diagonal terms of the interaction op- 




0.1 0.5 1 

sin(8/2) 



FIG. 5: (Color online) The decay of the edge Green's function 
corresponding to the inset 3-electron occupation numbers. As 
previously, the short red lines illu strate slop es —g = —3 and 
— 1, and Ge is evaluated at r = \/3(N — 1) + 1. 
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(a) (b) (c) 

FIG. 6: The pair-correlation function g{z\,z%) — 
pz{zi, Z2)/p(zi) with the first coordinate placed at the den- 
sity maximum of the negative y-axis zi = for (a) exact 
diagonalization, (b) 2-RDMFT, and (c) 1-RDMFT with P- 
0.675. Contours are separated by 0.05 l/ivl 2 and start from 
0.05 1/tvI 2 in (a) and (b) and from in (c). 



erator. The 1-RDMFT that only uses the diagonal Vijij 
and Vijji terms of the interaction matrix can not yield 
the correct behavior unless we have a very good density 
matrix functional. 

Recall that the backbone of the 1-RDMFT was the 
approximation of the pair-density. Figure [6] shows the 
pair-correlation functions corresponding again to exact 
diagonalization, 2-RDMFT, and P-0.675, where the lat- 
ter is reconstructed from the 1-RDM using Eq. Q. Al- 
though the 1-RDM reconstructed pair-correlation func- 
tion is reasonable vanishing at the position of the fixed 
electron, only the 2-RDMFT is able to produce the two- 
peak structure of the exact result with reduced density 
along the y-axis. 

Granted that the interaction energy functional in the 
2-RDMFT is exact, the results still do not coincide 
with the exact diagonalization results because the 3- 
representability conditions that ensure the physicality in 
this 3-electron system can not be taken into account with- 
out invoking the generalization of 2-RDMFT to include 
higher order RDMs. Figure [7]^a) shows the exact non- 
zero matrix elements of the 2-RDM Tny) m th° antisym- 
metric basis = — \ji))/2 and (b) calculated 
with the 2-RDMFT. The largest discrepancy between the 
results is the vanishing of the matrix elements involving 
states |(14)) or | (25)) for the exact result, while for ex- 
ample r^j ^ in the 2-RDMFT result. The matrix 
elements should vanish, as they are related to the ficti- 
tious many-body basis states that have double occupancy 
of orbital 4 or 2 (1 + 4 + 4 = 2 + 2 + 5 = 9, the total angu- 
lar momentum) . Consequently, the absolute values of the 
matrix elements also differ. However, the computed pair- 
correlation function and edge Green's function suggest 
that many of the physical quantities are not significantly 
affected by the absence of exact A^-representability, and 
that the lack of computer power might be the only real 
stumbling stone in the way of the 2-RDMFT. 



VI. CONCLUSIONS 

In summary, we have applied the one-body reduced 
density-matrix functional theory to small and large quan- 
tum Hall droplets at the spin-polarized strong magnetic 
field regime. The density-matrix power functional seems 
to work reasonably well at the strongly correlated v <C 1 
regime where the occupation numbers of the natural or- 
bitals are small. The newly applied method yields pre- 
viously inaccessible valuable information with large par- 
ticle numbers about the energetics and quantities that 
derive from the one-body reduced density matrix. The 
density-matrix power functional yields reasonable bulk 
densities with the power parameter in the range 0.65- 
0.7. However, the detailed properties of the edge are not 
produced accurately with this functional. Moreover, it 
is not known if a good functional for a specific quan- 
tum Hall state would work universally at different filling 
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FIG. 7: (Color online) Illustration of the 2- RDM for three electrons in the v = 1/3 state computed with (a) exact 

diagonalization and (b) 2-RDMFT. The same color bar applies to both figures while zero values are left white. The negative 
matrix elements are indicated by black edge. 



fractions. Nevertheless, the computationally expensive 2- 
body reduced density matrix method seems to facilitate 
the properties of the edge, though this should be verified 
with a larger electron number in future. 

Prospects of the 1-RDMFT in quantum Hall systems 
include generalizations to spin and multiple Landau lev- 
els. Studies with systems without edge (sphere) and non- 
trivial topology (torus) are also encouraged while new 
state of the art energy functionals are of course very wel- 
come. 
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